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A model of globally coupled phase oscillators under equilibrium (driven by Gaussian white noise) 
and nonequilibrium (driven by symmetric dichotomic fluctuations) is studied. For the equilibrium 
system, the mean-field state equation takes a simple form and the stability of its solution is examined 
in the full space of order parameters. For the nonequilbrium system, various asymptotic regimes 
are obtained in a closed analytical form. In a general case, the corresponding master equations are 
solved numerically. Moreover, the Monte-Carlo simulations of the coupled set of Langevin equations 
of motion is performed. The phase diagram of the nonequilibrium system is presented. For the long 
time limit, we have found four regimes. Three of them can be obtained from the mean-field theory. 
One of them, the oscillating regime, cannot be predicted by the mean-field method and has been 
detected in the Monte-Carlo numerical experiments. 
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PACS numbers: 05.40.-a, 05.60.-k 



INTRODUCTION 



A system of coupled oscillators has been treated as 
a model system of collective dynamics that exhibits a 
plenty of interesting properties such as equilibrium and 
nonequilibrium phase transitions, coherence, synchro- 
nization, segregation and clustering phenomena. It has 
been used to study active rotator systems , electric cir- 



charge-density waves 
I , planar XY spin mod- 



cuits, Josephson junction arrays 
H oscillating chemical reactions 

els ||, networks of complex biological systems such as 
nerve and heart cells ||. 

Such a system of N-coupled phase oscillators is deter- 
mined by a set of equations of motion in the form [Q 

JV 

Xi=u>i + f(xi) + y^ j KijG{xj,x i ) + r)i(t), 
i=i 

i=l,...,N, (1) 

where Xi denotes the phase of the ith oscillator and uji 
is its local frequency, i.e. its frequency in the absence of 
the interaction between the oscillators. The local force 
is represented by the function f(x) and G{x,y) includes 
the coupling effect between oscillators. The constants 
Kij are the coupling strengths and rji(t) characterizes 
fluctuations in the system. In the case of weak coupling, 
G{x, y) = G(x — y) and G is a periodic function of its ar- 
gument. The specific model G(x) — sin a; has been inten- 
sively studied and in the physical literature it is known 
as a Kuramoto model Q. If Kij are positive then the 
coupling is excitatory (meaning Xi tends to pull Xj to- 
ward its value). If Kij are negative then the coupling is 
inhibitory (it tends to increase the difference between Xi 



and Xj). Most of studies of the model focus on the global 
coupling (each oscillator interacts with all the other os- 
cillators), where all pairs are interacting with uniform 
strength, Kij = K/N . Then the mean-field treatment 
holds exactly when N — > oo. 

In the paper we study a special case of the model 
(Q) when the fluctuation term represents thermal- 
equilibrium and nonequilibrium fluctuations. The re- 
mainder of this paper is organized as follows. In the next 
section we analyze an equilibrium system. It is a model 
with thermal fluctuations being Gaussian white noise. 
In Sec. [II , we study a nonequilibrium system by adding 
the second fluctuation source, i.e., a zero- mean, exponen- 
tially correlated symmetric two-state Markov process. It 
can describe a case when local fre quenc ies u>i of the os- 
cillators fluctuate in time. In Sec. [II B , we present the 
mean-field numerical solutions of a corresponding mas- 
ter equation and discuss results of the Monte Carlo sim- 
ulations of Langevin equations. Finally, in Sec. [V we 
formulate the main conclusions. 



II. MEAN-FIELD EQUILIBRIUM SYSTEM 

In this section, we analyze a system of phase oscillators 
in contact with thermostat of temperature T, namely, 



= - sinxi + — sin(a;j - Xi) + Ti(t), 



N 



N 



3 = 1 



' i v. 



(2) 



where thermal-equilibrium fluctuations Ti (t) are modeled 
by zero-mean delta-correlated Gaussian white noise, 

<r,(t)> = 0, (r^r^s)) - 2T6ijS(t - s). (3) 

This model can represent a planar model with anisotropy 
or external field. More general models than (0) has 



2 



been analyzed. Nevertheless, we reconsider the simpli- 
fied model (||) because of two reasons. Firstly, the state 
equation of the system has a simple tractable form. Sec- 
ondly, a new aspect of the stability problem of states is 
presented. 

Let us rewrite the interaction term in the form [|| 



1 N 

where the averages 



s cosxi — csmxi, 



1 



iV 



s = — > sin a% 



1 



N 



C = — > COS Xj 



(4) 



(5) 



are order parameters for the system ([!]). In the ther- 
modynamical limit, N — > oo, for each oscillator Xi — x 
the mean-field Langevin equation is obtained from the 
system (Fy) and reads 



x = F(x,s,c) +T(t), 



(6) 



where the effective force F(x,s,c) — —V'(x,s,c) (the 
prime denotes a differentiation with respect to x) and 
the effective potential 

V(x, s, c) = — (1 + Kc) cosx — Kssinx. (7) 

Let us introduce a probability density 

P(x,t) = (6(x(t)-x)) (8) 

of the process (0), where x{t) is a solution of (^|) for a fixed 
realization of noise T(t) and (...) denotes an average over 
all realizations of T(t). This density is normalized on a 
real axis, 

oo 

P(x,t)dx = l (9) 

and obeys the Fokker-Planck equation 

= £-V\x, s, c)P{x, t) + T^P(x, t). (10) 

The reduced probability density P(x, t) defined by the 
relation 



P(x,t)= ^2 P{x + 2-Kn,t) 



(11) 



n— — oo 



satisfies the Fokker-Planck equation ( jiOj ) as well, is peri- 
odic 

P(x + 27m, t) — P(x, t) for any integer n (12) 
and normalized on one period, 

x +2tt 

j P{x, t) dx = 1 for any real xq. (13) 



The order parameters s and c are determined self- 
consistently from the set of two equations H , 



oo 

s = (sin a;) = J sin xP(x.t) dx 



sin xP(x,t) dx = g(s, c), (14) 



oo 

cHoob*) = / - xP(x,t)dx 



cos xP(x, t) dx = h{s, c), (15) 



where P(x, t) = P(x, s, c, t) and P(x, t) = P(x, s, c, t) de- 
pend on parameters s and c via the effective one-particle 
potential V(x, s, c) given by (0). 

Our concern is the behavior of the system in the limit 
of long time, t — > oo. The stationary state is a thermo- 
dynamic equilibrium state and the stationary solution 
P s t(x) of ( [10| ) is a Gibbs distribution, 



P st (x) = Nc -V(x,s,c)/T^ 

N- 1 = J e - y ^ s ^' T dx. 
o 



(16) 
(17) 



It is well known that in the equilibrium state the average 
angular velocity vanishes (the principle of detailed bal- 
ance holds) (x) — (see, e.g., Q). Then from (||) it 
follows that 



(sinx) = 



(18) 



in the stationary state and only a symmetric state is real- 
ized for which the effective potential reduces to the simple 
form 



V(x, s, c) = V(x, 0, c) = — (1 + Kc) cosx. 



(19) 



The form of this potential is the same as for a system 
of non-interacting oscillators, V(x) — — cosx. However, 
the amplitude A = 1 + Kc can change. If c > then 
the coherence effect occurs and the most probable state 
is the deterministic state x = 0. On the other hand, if 
1 + Kc < then the most probable state changes and the 
new state is x = ir. 

The order parameter c is determined by the equation 



c/n 



1 + Kc 
T 



= h 



1 + Kc 
f 



(20) 



where Iq{z) and I\{z) are the modified Bessel functions. 
This equation can possess one, two or three solutions (see 
Fig. [l]). If the coupling strength K < 1 then only one 
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FIG. 1: The order parameter c = (cosx) for the system of 
coupled phase oscillators in equilibrium as a function of the 
coupling strength K and temperature T. All data have been 
obtained as solutions of the implicit equation (pp|). The up- 
per plot shows the dependence of c on the couplin g K for 
selected temperatures. Even at T — the equation (BOJ) has 
three solutions c = (ei, 02,03) — (1, —1, —1/K) for K > 1 (for 
K < 1 there exists only one solution ci = 1). The lower plot 
shows the order parameter c as a function of temperature. As 
it could be expected, for large thermal fluctuations, stochas- 
tic forces overwhelm the potential and the coupling, and the 
solution tends to ci = as T — > 00. The stability analysis 
shows that stable are only solutions with c > 0. 



solution exists for any temperature T of the system (Fig. 
[l]a). For high temperature, T » 1, the upper branch 
ci tends to zero as c ~ T _1 (Fig. [yb). The opposite 
asymptotics, when T — ► 0, can be obtained as well. In 
this case the upper branch c\ — ► 1 for any K > 0. The 
lower branch ci — ► — 1 and the middle branch C3 — > — 1/K 
for if > 1 (Fig. |l^i). Now, let us study stability of the 
stationary solutions. The linear stability analysis should 
be performed on the full set of equations of motion for 
average values s and c, (O) and ([l5|). Multiplication of 
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FIG. 2: The plot shows a vector field of the right hand side 
of the dynamical system (§|)-(|27|) for K = 3.1 and T — 0.5. 
Three black dots are solutions of the mean-field problem, i.e. 
stationary points of (pi|)-(|2^). The upper solution ci is a sta- 
ble node, the middle one C3 is an unstable node and the lowest 
one C2 is a saddle point. Let us notice that stability analy- 
sis in one dimension (assuming that s = 0) would lead to a 
false conclusion that the lowest point C2 is stable. Solid lines 
are a result of the Monte Carlo simulation of 8000 particles 
with the initial condition set to the mean-field solutions ci 
and C3. In the (s, c)-space the system evolves along the clock- 
wise or anticlockwise "semicircles" (depending on the initial 
microstate) to the stable node (0, ci). 



@ by either sin x or cos x and integration over x gives 

s = -(1+Kc) <sc> +Ts + Ks < c 2 >, (21) 
c = (1+Kc) <s 2 > -Tc- Ks <sc>, (22) 

and < ... > stands for the averages of products of cos a; 
and/or sina; (e.g. < sc >=< sinscosa; >). To make 
the system (|2l|)-(|2^) closed, we should write equations 
of motion for the unknown statistical moments < sc >, 



< s 2 > and < 



>. New, higher-order moments will 



occur and in this way we obtain a hierarchy of infinite 
number of differential equations for moments, which is 
difficult to handle. Therefore we proceed in another way. 
Let us notice that for < s 2 >=< sin 2 x > one may write 
< s 2 >= 1- < cos 2 a^ >= 1- < c 2 >. Additionally, 
one can introduce deviations from the mean values and 
write < c 2 >= c 2 + < (5c) 2 > as well as < sc >= s c+ < 
5s5c >. As a result, one obtains 

s = -(c — K < (Sc) 2 



> —T) s 
-(1 + Kc) < SsSc >, 
(1 + Kc) (1 - c 2 - < (Sc) 2 >) - Tc 

—K s < SsSc > 



Ks 2 c 



(23) 



(24) 
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From ( |i8| ) we know that s — is a stationary solution 
of the above equations. In order to obtain this solution 
s = from (|2^), the correlator < SsSc > should vanish 
in the stationary limit, i.e., < SsSc > — > as t — > oo. 
Insertion of s = into the second equation ( p4| ) with c = 
yields stationary solutions for c. They are determined 
by the equation 

[l + Kc)(l-c 2 - <{5c) 2 >)-Tc = (25) 

These solutions depend on the unknown variance < 
(6c) 2 >. In the low temperature limit T — > 0, the 
variance < (<5c) 2 >— > and we recover the solutions 
ci = 1, C2 = — 1 and C3 = —1/K. In this case, the linear 
stability analysis of (^)-(|24|) shows that the stationary 
point (0, ci) is a stable node, the point (0, C2) is a saddle 
and the solution (0, C3) is an unstable node. For T > 0, 
the stability of solutions remains unchanged. Indeed, in 
our simulations we have confirmed this statement. We 
have also analyzed an auxiliary dynamical system de- 
fined by a set of two differential equations, namely (cf. 
(0) and (|iT 



s = -s + g(s,c), 
c = —c + h(s, c). 



(26) 
(27) 



The stationary solution of this system is the same as the 
equilibrium state of the system (g). In Fig. ||we present a 
vector field generated by the dynamical system (|26|)-(|2~7|) 
and its three stationary points (sj,Cj),i = 1,2,3. One 
can infer that the upper point (0, c\) is a stable node, the 
lower point (0, C2) is a saddle and the middle point (0, C3) 
is an unstable node (the same as for (p3f)-(p4|)). We have 
also found unexpectedly that the trajectory of the system 
(^)-(^) is the same as that obtained from simulations of 
the set of Langevin equations (|J), see Fig. |^. It allows us 
to formulate the conjecture that the hierarchy of infinite 
number of equations for moments of the set (sin a;, cos x)is 
equivalent to (p6|)-(p7|). Unfortunately, we cannot prove 
it rigorously. 



Markovian stochastic processes [[L0| , 

&(*) = a > 0, 

P(—a — ► a) = P(a — > —a) = fx, 



(29) 



where P(— a — > a) is a probability per unit time of the 
jump from the state —a to the state a. This process is of 
zero average, (£t(i)) = 0, and exponentially correlated, 



Um(s))=a 2 S ij e 



-|i-a|A 



(30) 



where t = l/2/i is correlation time of the process £i(t). 
So, it is characterized by two parameters: its amplitude 
a (or equivalently the variance < £ 2 (t) >= a 2 ) and the 
correlation time r. 

The mean-field Langevin equation takes the form 



x = -V'(x,8,c) + T(t) + £(t) 

and the corresponding master equations read jflj 

dP+(x, t) d rT „. . . „ . 
g t - - ^ [V^x, s, c) - a] P+(z, t) 

+T^P+(x, t) - fiP + (x, t) + tiP_(x, t) 
dP_(x,t) d 



(31) 



dt 
d 2 



dx 



[V'(x,s,c)+a]P_(x,t) 



+T-^P- (x, t) + M P+ (x,t)- (z, «) 



(32) 



(33) 



where the probability densities 
P+(x,t) = P(x,a,t), P-{x,t) 



P{x,-a,t). (34) 



depend on the order parameters s and c, which in turn de- 
pend self-consistently on the marginal density P(x, t) — 
P+OM) + P-(x,i) via the relations @-(||). Eqs (||)- 
( p3| ) cannot be solved analytically, even in the stationary 
state. However, in some limiting cases, stationary so- 
lutions of them are known, e.g., if the correlation time 
t — > 00 (the adiabatic limit) or if temperature of the 
system is zero, T = 0. 



III. MEAN-FIELD NONEQUILIBRIUM 
SYSTEM 

Nonequilibrium systems can be modeled by including 
a term which describes non-thermal and nonequilibrium 
fluctuations, noise and perturbations. There are many 
possibilities to do this but here we consider a slight mod- 
ification of the previous model, namely, 

K N 

±i = - sinxi + — 2J sin(£j - x t ) + Ti(t) + ^(t), 

i = l,...,N, (28) 

The random functions £i(t) represent nonequilibrium 
fluctuations and are modeled by a symmetric dichotomic 



A. Analytical results 

From the ratchet theory we know that the stationary 
average angular velocity is zero, (v) = (x) = 0, because 
the potential (Q) is symmetric and fluctuations (^) are 
symmetric |l|. Therefore 

s = (sin a;) = (35) 

and V(x, s, c) takes the same form as in the previous case 
(|l9|). In the adiabatic limit, the equation determining a 
stationary state is 

c = — J cos x \p-\-(x, c) +p-(x, c)] dx, (36) 
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where the stationary probability densities 
Pi(x,c), (i = +, -) read 



Pi(x,c) = 



C/ l (x,c)" : } 2 Vr 1 (y, c )dy 

X 

2-K X + 27T 

J ' Ui(x,c) J U~ 1 (y,c)dydx 

x 



(37) 



and 



U±(x,c)=e v ^°< c V T e ±ax / T . 



(38) 



In the second limit, i.e. when temperature of the sys- 
tem is zero, T = 0, the stationary state is determined by 
the equation 



/ cos a; J D- 1 (x,c)e-*( ;E ^) dx 

= n(c) 

/ D- 1 {x,c)e-^( xc '> dx 
n(c) 

where the thermodynamic potential 

¥(s,c) = J D- 1 (y,c)V'(y,0,c)dy 
o 

and the effective diffusion function 

D(x,c) = t [a 2 - V'(x,0,c) 2 ] . 



(39) 



(40) 



(41) 



The integration interval fi(c) = [0, 27r] iff D(x, c) > 0. 
If in some intervals the function D(x, c) is negative then 
f2(c) = [xi,X2], where x\ and X2 are suitable roots of 
the equation D(x,c) = and in the interval [xi,^] the 
diffusion function is positive. 

The limiting case T = and r — > oo is analytically 
tractable. From the master equations it follows that in 
this case the stationary state is determined by the equa- 
tion 



[a 2 - V'(x, 0, c) 2 ] P(x) = const. 



(42) 



In the diffusive regime, when dichotomic noise activates 
both forward and backward transitions over barriers of 
the effective potential, the solution of (E4) is 



P(x) = const./ [a 2 - V'(x, 0, cf 



(43) 



and the only solution of the state equation (|15| ) is c = 
0. In the non-diffusive regime, when dichotomic noise 
cannot activate neither forward nor backward transitions 
over barriers of the effective potential, the normalized 
solution of (^3) has the form 

P(x) = ^[S{x-x 1 ) + S(x-x 2 )} (44) 

where X\ and x 2 are solutions of the equation 

a 2 - V'(x,0,c) 2 = 0. (45) 



If 1 + K c > then the state equation is determined by 

c = cos [arcsin(a/ (1 + Kc))] (46) 

This equation can possess two positive roots, c x > c 2 > 0. 
The solution c\ is stable while c 2 is unstable. If 1 + Kc < 
then 



c = — cos [arcsin(a/(l + Kc))] 



(47) 



This equation can possess two negative roots which are 
unstable. It is depicted in Fig. |[ 



B. Numerical methods 

In a general case the mean-field problem reduces to the 
set of non-linear master equations ([32]) and ( |33| ) which 
has to be solved. Apart of a few previously considered 
special cases which can be treated analytically, only nu- 
merical methods are applicable. We have approached the 
numerical problem of solving ( |32] ) and (|33) as follows. 
Conditions of self-consistency have been considered ( |l4| ) 
and (15) as the non-linear minimization problem in two 
dimensions on the bounded domain — 1 < s < 1 and 
— 1 < c < 1. It has been handled in a standard way, 
making use of numerical libraries. However, each eval- 
uation of functions g and h for given (s, c) requires the 
knowledge of stationary solution of the system ( ^2|) , (|3^ ) 
with fixed s and c. This is, in turn, a linear boundary 
value problem which can be easily solved with the help 
of finite element method (FEM). In the case T = 0, the 
stationary distribution P(x) is given be quadratures but 
it is very difficult to handle it analytically. Moreover, 
we found it practically easier to obtain a FEM solution 
for small enough T then to estimate, divergent in some 
cases, triple integrals. 

Additionally, in order to verify mean-field results we 
have performed Monte-Carlo simulations of the Langevin 
equations (j28|). This, independent method enabled not 
only verification of numerical results but also applicabil- 
ity of the mean-field approach. Because the Monte-Carlo 
simulation follows the evolution of microscopic state of 
the system it can be considered as a numerical exper- 
iment, in contrary to the mean-field approach which is 
only an approximation. In general, Monte-Carlo simu- 
lations of globally interacting iV-particles require ~ N 2 
operations per time step. The special form of the in- 



teraction term ~ sin(xj — Xi), leads to relations (^]) and 
(|). In the course of simulation the average values s and 
c need to be evaluated only once per a simulation step, 
what reduces the number of operations per a time step 
to ~ N. 



C. General case: numerical analysis 

All numerical mean-field results have been obtained 
in the stationary regime. First we study the zero- 
temperature case, T — 0. The natural characteristics 
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x=2.0 



C O^ffss 




FIG. 3: Solutions of the stationary mean-field problem 
([bi|),(|l5|), ( |3^ ) and (p3|). The temperature is zero and the 
amplitude of dichotomic fluctuations is a = 2.8. As in the 
equilibrium case, negative solutions are unstable. Stable, and 
observed in Monte-Carlo simulations solutions lie in the up- 
per plane c > 0. Grey regions depict places where the locking 
condition is met i.e. maximum of the effective force over- 
whelms the amplitude of fluctuations. One can notice that 
for small values of r, the system starts to behave as an equi- 
librium one. In the case r — * oo the asymptotic analytical 
solution is shown. For finite r numerical results are depicted. 
In this case the system exhibits onset of hysteresis in c(K). 



of the stationary state are statistical moments, in partic- 
ular the first two moments < x > and < x 2 >. These 
moments are not good characteristics in the case consid- 
ered. If the system is spatially periodic, then for any 
spatially periodic function A{x) = A(x + 2tt) we can cal- 
culate its mean value exploiting either the probability 
density P(x, t) or the reduced probability distribution 
P{x, t) because then the equality 



2tt 



(A(x)) 



fective potential (19)). These two regimes, marked by 
white and gray regions in Fig. |[ are separated by two 
critical lines: Kc + 1 = a for positive values of c and 
Kc + 1 = —a for negative values of c. For negative c, 
the dependence of c upon K is qualitatively the same 
as for the equilibrium system (Fig. [j]). These solutions 
are unstable and therefore will not be considered. Now, 
let us discuss the positive solutions c > 0. They depend 
strongly on the correlation time r of dichotomic fluctu- 
ations. For short correlation time, the order parameter 
c monotonically increases with the growing coupling K. 
For longer correlation time, new effects arise: the depen- 
dence is discontinuous and hysteretic. In some domain 
there are three solutions ci > C2 > C3. The solutions 
ci and C3 are stable while C2 is unstable. The hystere- 
sis is bigger and bigger when r increases. The jumping 
point K\ from the lower to the upper branch tends to 
infinity and the jumping point K2 from the upper to the 
lower branch tends to a constant value determined by 
eq. (|46|). The upper branch of solutions C\{K) — > 1 and 
the lower branch c^{K) — > when r — > 00. For r = 00, 
the solutions split into two branches of three solutions, 
namely, one C3 = and two other determined by fl46|). 
The stationary mean-field solutions have been verified 
by the Monte-Carlo simulations. The comparison is pre- 
sented in Fig. 0. Simulations show that the implicit as- 



0.75 

C 

0.5 



0.25^ 



K=2.94 




I A-*. El 



A(x)P(x, t) dx = A(x)P(x,t) dx (48) 



3 

K 



holds. It is not a case for non-periodic functions and then 
there is a problem which the distribution should be used 
for calculationg the average value. Therefore we consider 
periodic functions. Here, two natural order parameters 
s =< sinx > and c =< cosx > occur which characterize 
the probability distribution in the same way as < x > 
and < x 2 >. Indeed, the function sin a; is odd like the 
function x and the function cos x is even like the function 
x 2 . Because < sins >= 0, below we analyze < cosx >. 
In Fig. H we show the dependence of the order param- 
eter c = (cosx) on the coupling strength K. One can 
distinguish two main regimes: the diffusive (dichotomic 
noise activates transitions over barriers of the effective 
potential (|l9|)) and non-diffusive or locked (dichotomic 
noise cannot activate transitions over barriers of the ef- 



FIG. 4: The same case as in Fig. ^: the comparison of mean- 
field results and Monte-Carlo simulations. For r = 0.25 and 
r = 0.5 there is a perfect agreement of the Monte Carlo and 
mean-field methods. However, for r = 1 and r = 2 temporal 
oscillations of density of particles appear in the Monte-Carlo 
method. Thus the order parameter also performs temporal 
oscillations. An averaged value of those oscillations differs 
from those coming from the mean-field solution. The stan- 
dard deviation of c (averaged in time) is shown in the lower 
insert. One can notice that if oscillations disappear (Sc — 0) 
then the simulated values of c agree very well with mean-field 
predictions. 

sumption of time-independent stationarity of the system 
(when t — > 00) is restricted to some values of parameters 
of the model. Indeed, if the time-independent station- 
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ary state of the system exists then the mean-field solu- 
tions agree with simulations. In particular, for r = 0.5 
the hysteresis is observed (see point K — 2.94 in Fig. 
|] and Fig. 0). However, for longer correlation time t, 
temporarily oscillating steady-states exist for which the 
probability distribution P(x, t) is periodic in time. In 
consequence the order parameters s — s(t) and c = c(i) 
are time-periodic and in the limit of long time, the time- 
dependent steady-states appear. This is the case when 
the mean-field predictions fail, e.g. the hysteresis is not 
realizable. In Fig |^ we depicted this phenomenon for 
t = 1,2. We have noticed only monotonic dependence 
of c upon K (if c is a periodic function of time, its 
time average is taken). The quantity which can charac- 
terize the time-independent stationarity/time-dependent 
stationarity (i.e. oscillations) of the long time state is the 
time-averaged standard deviation (5c) = (c 2 )t — (c)f of 
the order parameter. We have observed that if (5c — 
then the mean-field solutions are correct. Otherwise, 
they are incorrect. It is shown in the lower insert of 
Fig. g 



1.5 



T 1 



0.5 



0, 




Unlocked 



1 



FIG. 5: The phase diagram of the system with a — 2.8 and 
T — 0. Five various regimes are distinguished: unlocked and 
oscillating, mean-field hysteretic, locked and bistable. The os- 
cillating regime has been verified by Monte-Carlo simulations. 
All other data come from the stationary mean-field problem. 
The empirical formula r = 2/K surprisingly well fits the left 
boundary of the oscillating region. 

In Fig. U we present the phase diagram on the (K, r) 
plane for a fixed amplitude a = 2.8 of dichotomic fluctua- 
tions. In the case of non-interacting oscillators, this value 
of a corresponds to the diffusive regime. Roughly speak- 
ing, there are two regions: diffusive when a > 1 + Kc(t) 
(i.e. dichotomic noise can induce transitions over bar- 



riers of the effective potential (19)) and non-diffusive 
when a < 1 + Kc(t) (i.e. dichotomic noise cannot in- 
duce transitions over barriers of the effective potential 
( p"9| ) ) . The diffusive region is divided into two parts which 
we call the unlocked regime (where the mean-field solu- 



tions are correct) and the oscillating regime (where the 
mean- field solutions fail). In the unlocked regime, there is 
one and only one time-independent stationary state and 
there is only one stationary value of the order parame- 
ter c =< cos a; > which is always stable. In this regime, 
the reduced stationary probability density P s t(x) ^ 
for any x. It means that with non-zero probability the 
phases of the oscillators can take any value of x and os- 
cillators are not synchronized. In the oscillating regime, 
the only stationary state is temporarily oscillating state 
for which lim^oo P(x, t) is time-periodic and the order 
parameter c = c(t) is time-periodic. From previously 
discussed results it follows that this regime is bounded 
from the right, i.e. if K > Kq then this regime disap- 
pears. The critical value Kq can be determined by Eq. 
( pitf ) from the condition that it possesses the double root 
c\ = C2- The oscillating regime is presented in Fig. [(], 
where we show evolution of two distributions, the full 
density P(x,t) normalized on the interval (— oo, oo) and 
the reduced density P(x, t) normalized on one period. In 
the latter case, the density oscillates between the distri- 
bution of one maximum around tt (it corresponds to the 
maximum of the local potential — cos x) and the distri- 
bution of two maxima around and 27r (it corresponds 
to the minima of the local potential). 



P(X,t) 




671 o 



FIG. 6: Monte Carlo simulations of the system for a = 2.8, 
T = 1.0, K = 2.94 and T = 0. None of the mean-field predic- 
tion is realized. The only stable solution is a stationary one. 
In the upper plot evolution of the probability density reduced 
to i £ (0, 2tt) is shown. In lower plot the full distribution 
is presented. A starting value was a uniform distribution on 
x e (0,2vr). 

In turn, the non-diffusive region is divided into two 
other parts which we call the locked and hysteretic 
regimes. In the locked regime, only one steady-state 
solution exists. In this regime, there are intervals of 
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x for which the reduced stationary probability density 
Pst(x) = and the phases of the oscillators are locked 
in these intervals. It is an effect of interaction and cor- 
responds to the synchronization of oscillators (let us re- 
member that in the case on non-interacting oscillators, 
the diffusive regime is realized in which the phases can 
take arbitrary values). The synchronization is stronger if 
the support of P st (x) is smaller. In this regime, there is 
only one mean-field value of c =< cos x > which is always 
stable. The so-called MF hysteretic regime is defined in 
the following way: There are three mean-field stationary 
values ci > C2 > C3 > of the order parameter. The so- 
lutions ci and C3 are stable while ci is unstable. However, 
in this regime only one mean-field solution ci is realized, 
which lies on the upper branch of the mean-field hystere- 
sis, cf. the case r = 2 for K > 4 in Fig. [|. There is also 
a regime of bistability. As in the previous case, there are 
three mean-field stationary values c\ > ci > C3 > of the 
order parameter. But now, two stable solutions C\ and 
C3 can be realized what is demonstrated in Fig. ^. The 
upper state C\ > C3 corresponds to the locked regime, 
a < 1 + Kci(t) and the lower state C3 corresponds to the 
unlocked regime, a > 1 + Kc^{t), cf. Fig. ||, the case 
K = 2.94 and r = 0.5. It is also instructive to see how 




of all "particles" and realizations of noises are different, 
it determines evolution of P(x,t). For illustrating ani- 
mations of the time evolution we refer to our webpage 
& 

The influence of temperature is depicted in Fig. M (only 
the mean-field case is shown). On the basis of these 
results, one may conclude that the increase of thermal 
fluctuations acts like the decrease of correlation time r 
of nonequilibrium fluctuations. The hysteretic region in 
K is reduced as temperature grows. In particular in Fig. 
|we see that for T = 1.5 the mean- field problem has got 
only a single solution. 




K 

FIG. 8: The order parameter c versus coupling strength 
K for selected values of temperature T. The increase of T 
decreases the region of hysteresis. Remaining parameters are 
r = 2.0 and a = 2.8. 



FIG. 7: Monte Carlo simulations of the system for a — 2.8, 
t = 0.5, = 2.94 and T = 0. The starting point was 8000 
particles distributed uniformly on a; £ (0, 2tt). The only differ- 
ence between left and right scenarios is the microscopic state: 
individual particles were chosen differently (all macroscopic 
parameters are the same). The left scenario leads to diffusive 
state i.e. c < (2.8 — 1)/K while the right one leads to locked 
one c > (2.8 — 1)/K- The shape of a stationary, mean-field 
distribution is shown for t — > 00. 

the probability distributions P(x,t) or P(x, t) evolve in 
time approaching the long time limit. In Fig. |7j, the evo- 
lution of the density P(x, t) is shown for the values of pa- 
rameters chosen from the bistability regime of the phase 
diagram, i.e., when two stable stationary solutions exist. 
One can observe that in dependence of the microscopic 
initial conditions the system evolve either to the diffusive 
stationary state or to the non-diffusive locked stationary 
state. In two cases, the macroscopic state, i.e., the ini- 
tial probability density of oscillators is the same uniform 
distribution. The microscopic state, i.e., initial positions 



IV. SUMMARY 

In this paper we have investigated the equilibrium and 
nonequilibrium system of coupled phase oscillators. In 
fact, it can be any abstract model of interacting parti- 
cles in spatially periodic structure with a periodic global 
interaction (e.g. interacting Brownian motors jl4], 
The equilibrium system defined by eq. (^|) is a special 
case of models considered in the literature. Nevertheless, 
to the best of our knowledge, the state equation fl2C| ) 
has not been presented. We pay attention to the sub- 
tle stability problem which sometimes is treated superfi- 
cially |l5| . Properties of the nonequilibrium system fl2S| ) 
are naturally much more interesting. The phase diagram 
consists of five parts and cannot be fully obtained from 
the mean-field approach. The non-mean-field regime is 
the oscillating regime, which has been detected by use 
of the Monte-Carlo simulations and by analyzing fluc- 
tuations of the order parameter c = (cos a;). The next 
interesting finding is that although the non-interacting 
system is in the diffusive regime, the interaction can 
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move the system to the non-diffusive regime and then 
"particles" are confined in valleys of the potential (of 
course it is exact when temperature T = 0). It means 
that effectively, for the one-particle dynamics, the barrier 
height 2(1 + Kc) of the local potential is magnified and 
nonequilibrium fluctuations of amplitude a are not able 
to induce transitions over barrier. 

All the results so far refer to the simple reflection- 
symmetric local potential — cos:c. If we add the higher 
order harmonics, e.g. cos 2x, the potential is still sym- 
metric. However, behavior of the system can then be 
radically different because the second order parameter 



s = (sinx) =/= 0. New phenomena such as the symmetry 
breaking, phase transitions and noise-induced transport 
can occur in the system. The paper on this subject will 
be published elsewhere. 
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